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Working in the crystal-momentum representation, we calculate the optical conductivity of non- 
centrosymmetric insulating crystals at first order in the wave vector of light. The time-even part of 
this tensor describes natural optical activity and the time-odd part describes nonreciprocal effects 
such as gyrotropic birefringence. The time-odd part can be uniquely decomposed into magnetoelec- 
triclike and purely quadrupolar contributions. The magnetoelectriclike component reduces in the 
static limit to the traceless part of the frozen-ion static magnetoelectric polarizability while at finite 
frequencies it acquires some quadrupolar character in order to remain translationally invariant. The 
expression for the orbital contribution to the conductivity at transparent frequencies is validated by 
comparing numerical tight-binding calculations for finite and periodic samples. 
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I. INTRODUCTION 

Electric and magnetic effects are closely coupled in 
magnetoelectric (ME) materials. These are insulators 
with broken spatial-inversion (V) and time-reversal (7") 
symmetries, in which an applied electric field £ induces a 
first-order magnetization M, and conversely a magnetic 
field B induces a first-order electric polarization P. This 
cross response is described in the static limit by a single 
magnetoelectric polarizability tensor 
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where the equality follows from changing the order of the 
mixed derivatives of the free energy. 

The ME effect has been intensively studied in recent 
years. While the focus has been mostly on the static 
response, ME effects in the optical range have also been 
observed. 1 For oscillating fields the thermodynamic argu- 
ment leading to the second equality in Eq. (1) does not 
hold because the system is not in equilibrium, and two 
separate frequency-dependent polarizabilities are needed 
to describe the dynamical ME coupling 
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It was recognized already in the 1960s that the coupling, 
Eq. (2), leads to new optical effects in ME media, such 
as gyrotropic birefringence. 2 Since the lattice-mediated 
response is frozen out at optical frequencies, the purely 
electronic contribution can be isolated. The first success- 
ful measurements, on Cr2C>3, found that the strength of 
the optical ME coupling is comparable to that of the 
static one. 3 

The phenomenology of optical ME effects has been 
studied in detail in the literature, starting with the work 
of Hornreich and Shtrikman on gyrotropic birefringence. 4 
These authors showed that this effect is a consequence of 
spatial dispersion, appearing at first order in the expan- 
sion of the effective optical conductivity tensor (defined 



by Eq. (6) below) in powers of the wave vector q of light 
a ab (q_,u}) = a^ b '(u)) + a abc (oj)q c H (3) 

It is well known that the phenomenon of natural opti- 
cal activity is also a manifestation of spatial dispersion. 5 
While natural optical activity is associated with the T- 
even part of (J abc (uj), optical ME effects arise from the 
7~-odd part, which can be nonzero only in magnetically 
ordered systems, where T symmetry is spontaneously 
broken. A careful consideration of all response tensors 
which contribute to the conductivity at linear order in q 
shows that these include, in addition to the dynamic ME 
polarizabilities, Eq. (2), the electric-quadrupole response 
of the medium. 

Regarding the microscopic theories needed for quanti- 
tative calculations, there are well-established molecular 
theories of spatial dispersion, 6,7 but the corresponding 
theory for crystals is not equally developed. A band the- 
ory of natural optical activity was put forth by Natori 8 
but has not been used in first-principles calculations. To 
our knowledge, only one group has reported calculations 
of natural optical activity in solids at optical wavelengths, 
based on a somewhat different formulation. 9,10 As for the 
optical ME effects, quantitative estimates of their mag- 
nitude have so far relied on cluster models to mimic the 
crystalline environment. 11,12 

In this work, we develop a formalism for calculating 
spatial-dispersion effects in the framework of band the- 
ory. One difference with respect to previous works is 
that we give a unified treatment of both T-even and 
T-odd parts of this tensor. More importantly, we ex- 
press the transition matrix elements in the crystal mo- 
mentum representation. 13 This choice has both practical 
and formal advantages. The practical advantage is that 
it leads to expressions which can be easily implemented 
using localized Wannier orbitals. On the theoretical side, 
the crystal- momentum representation is the language in 
which the modern theories of electric polarization, 14,15 
orbital magnetization, 16-19 and orbital magnetoelectric 
response 20,21 are formulated. As we shall see, our ex- 
pression for the orbital contribution to the T-odd part of 
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Cabc(w) generalizes to finite frequencies the traceless part 
of the orbital ME polarizability formula of Refs. 20 and 
21. 

The manuscript is organized as follows. In Sec. II we 
give a self-contained account of the phenomenology of 
spatial-dispersion optics. The effective conductivity is 
defined and related to the magnetoelectric and quadrupo- 
lar polarizabilities. We then reformulate the phenomcno- 
logical relations, originally obtained for finite systems, in 
terms of translationally invariant quantities which remain 
well defined in the thermodynamic limit. The main re- 
sults of the paper are contained in Sec. Ill, where we 
obtain a microscopic expression for the cr abc (uj) in peri- 
odic insulators. We then consider the u> — > limit of that 
expression and discuss its relation to the theory of static 
ME response. In Sec. IV we implement the bulk a a bc{u) 
expression for a tight-binding model and compare the re- 
sults with calculations on finite samples cut from the bulk 
crystal. We conclude in Sec. V with a brief summary and 
outlook. 



II. PHENOMENOLOGY OF SPATIAL 
DISPERSION 

In this section we discuss spatial dispersion from a 
phenomenological perspective. Besides introducing basic 
definitions and setting the notation, the main purpose 
here is to arrive at Eqs. (25)-(29) relating the spatially 
dispersive optical conductivity to translationally invari- 
ant renormalized multipole polarizabilities. Those rela- 
tions will allow us to identify the magnetoelectriclike and 
purely quadrupolar parts of the optical response of crys- 
tals, to be calculated in Sec. III. 



A. Effective conductivity tensor 

Consider a crystal with broken V and possibly broken 
T symmetries. We are mainly interested in materials 
where those symmetries are broken spontaneously, rather 
than by static electric and magnetic fields, and wish to 
study their current response J(q, w) to an electromag- 
netic plane wave 
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B(r,t) = -[q X £(q,cj)]e J ( q - r ^. 
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Because the oscillating electric and magnetic fields S 
and e and B are interdependent, the linear (in the field 
strengths) response can be described by a single effective 
conductivity tensor 4 ' 22 



J a (q,w) = cr a6 (q,w)£ b (q,w). 



(6) 



related (in Gaussian cgs units) by 

4:7rz 

£afc(q,w) = S ab H cr ah (q,w). 

UJ 



(7) 



The leading term in the expansion of a ab (q, uj) in pow- 
ers of q, Eq. (3), is the optical conductivity in the electric- 
dipole approximation. We shall focus on the next term 
in the expansion, cr abci which is chiefly responsible for 
spatial dispersion. Because spatial inversion takes q 
into q, the tensor <r abc (uj) necessarily vanishes in cen- 
trosymmetric systems. Its symmetric (a^ bc ) and anti- 
symmetric (er^, c ) parts under the interchange of the first 
two indices are, respectively, odd and even under T. 23 
The T-even piece describes natural optical activity, and 
the T-odd piece describes non-reciprocal optical effects. 
These include, in addition to gyrotropic birefringence, di- 
rectional dichroism 1 and magnetochiral effects in chiral 
ferromagnets. 24 

Unlike the spontaneous magneto-optical effects coming 
from the T-odd part of a^j (magnetic circular dichroism 
and birefringence), which require ferromagnetic or ferri- 
magnetic order, gyrotropic birefringence can also occur 
in antifcrromagnets such as Cr2 03. This is a well-known 
magnetoelectric material, and indeed the physical basis 
for spatial dispersion rests in part on the magnetoelectric 
effect. 



B. Multipole theory for finite systems 

The connection between spatial dispersion and the 
magnetoelectric effect can be readily established by ex- 
pressing J(q, ui) in terms of the multipole moments of the 
charge and current distributions. We begin by taking the 
spatial Fourier transform of the current density, 



J(q,t) = l / ,/r, ' q ' l 'J(r. /) 



(8) 



Alternatively, one may choose to work with the dielec- 
tric function e &(q, w). 5,22 To first order in q the two are 



and expanding in powers of q, 

J(q,t)=jW(t)+J«(q,t) + 0( g 2 ). (9) 

Standard multipole-expansion manipulations 22 involving 
the continuity equation and integrations by parts show 
that J ( a °\t) = d t P a {t) and 

jW(q,t) - - l -fd t Q ab (t) + ie abc cq b M c (t), (10) 

where e abc is the antisymmetric tensor of rank three and 
P, Q, and M are the electric dipole, electric quadrupole, 
and magnetic dipole moments of the sample divided by 
its volume 

P a (t) = ^J dvr a p{t,v), (11) 
Qab(t) = ^ J drr a r hP (t,r), (12) 
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M a (t) = ^y€ abc J drr b J c (t,r). (13) 7afcc = X ^c + xg ac X q ac b X\ ca 



Fourier transforming in time we arrive at 

J a (q,w) = -iuP a {u)-^qbQab{u)+i£abcCqbM c {Lu)+0{q 2 ). 

(14) 

The current induced by the monochromatic wave, 
Eqs. (4) and (5), can now be calculated from the oscil- 
lating induced moments, which are the real parts of the 
following expressions: 6,7 



Pa = X c ab £b + 2^ bc V c £ b + • • • + XafBb + • • • , (15) 



Qab 



X q abc£c + 



Xab tb 



(16) 



(17) 



where the fields and their gradients are evaluated at the 
location of the sample. x° is the electric polarizability per 
unit volume, and quantum-mechanical expressions for 
the remaining response tensors are listed in Appendix A. 
X cm and x mc ar e the dynamic ME polarizabilities intro- 
duced in Eq. (2); they involve matrix elements of the 
electric-dipole (El) and magnctic-dipole (Ml) operators, 
and for this reason are known as the El. Ml terms. x q 
and x q are the E1.E2 terms, as they mix electric-dipole 
and electric-quadrupolc transitions. 

In Eqs. (15)~(17) only those terms which contribute 
to the effective conductivity up to first order in q were 
kept. Combining Eqs. (14)-(17) with Eqs. (4) and (5) 
and comparing with Eqs. (6) and (3) we find, upon col- 
lecting terms linear in q, 



= habc + Xbac\ > 

(23) 



be 
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Xbac 



X q acb ~j~ X q bca j_ lp„r...q ...q 1 
^ - TjrtC [Xabc - Xbacl 



(24) 

In each of these equations the second equality, denoted by 
the symbol =, only holds at nonabsorbing frequencies, for 
which xlTHxZ)* and X q abcHxl ab Y (see Appendix A). 
In this lossless regime a a bc becomes anti-Hermitian in the 
first two indices. 

The above multipole formulation leads to a practical 
scheme for calculating spatial dispersion effects, by com- 
puting the polarizabilities x 0K \ X mc > X q > an d X q from 
Eqs. (A1)-(A4), and assembling them in Eq. (18). This 
approach can be used for molecules and other finite sys- 
tems but not for bulk crystals, because the quantum- 
mechanical expressions in Appendix A become ill-defined 
under periodic boundary conditions. 

The problem can be traced back to the integrations by 
parts carried out around Eq. (10), where the boundary 
terms were discarded. Such procedure is allowed for finite 
systems, as the boundary can always be placed outside 
the sample. It cannot, however, be rigorously justified for 
periodic crystals with delocalized electrons. This is a sub- 
tle but by now well-understood problem. For example, 
the macroscopic electric polarization and orbital magne- 
tization of crystals cannot be calculated under periodic 
boundary conditions as the first moments of the charge 
and orbital current distributions in one crystalline cell 
because the result depends on the choice of cell. 15 The 
correct band-theory expressions for P and orbital M have 
been derived in Ref. 14 and Refs. 16-19, respectively. 



Oabc = ic(Xa*d e dbc + 



cdXdb) + -^(Xlbc - X q acb) 



(18) 



Spatial dispersion is thus governed by the magnetoelec- 
tric and quadrupolar responses of the medium. 4 The need 
to include the quadrupolar terms in order to properly 
describe the optical activity of oriented molecules and 
uniaxial crystals was emphasized in Ref. 25. 

Dividing Eq. (18) into symmetric (magnetic) and an- 
tisymmetric (natural) parts under a «-» b yields 



a abc = ic [tbcdOtad + CacdCXbd) + ^Jabc, 



a tbc = * c ( e bcd/3ad - CacdPbd) + ^^abc, 



where we have defined 



a ab = Xab ^ ba =Rexg 



cm 
b ) 



which reduces to Eq. (1) in the static limit, and 



y em ...me 

Pab= ^ 2 Xba =XmXa 



cm 
b ) 



(19) 



(20) 



(21) 



(22) 



C. Translationally invariant polarizabilities 

Already for finite systems the description based on 
Eqs. (15)-(17) is highly redundant, as the individual po- 
larizabilities are origin dependent. 6 ' 7 The combination of 
polarizabilities on the right-hand side of Eqs. (19) and 
(20) is of course translationally invariant (the conductiv- 
ity is a physical observable) but we shall go one step fur- 
ther and redefine the polarizability tensors so that they 
become individually origin independent, and hence well 
defined for periodic crystals. 

To begin, we note that the trace of a drops out from 
Eq. (19), leaving eight magnetoelectric quantities. These 
fully specify a^ bc in the static limit while at finite fre- 
quencies the quadrupolar tensor j a b c — jbac contributes 
18 additional quantities. This brings the total number 
to 26, while a^ bc itself, being symmetric in the first two 
indices, only contains 18 independent quantities. The 
source of this discrepancy lies in the origin-dependence 
of the tensors a and 7, and it can be removed by suitably 
redefining them. To that end we note that any third-rank 
tensor a^ bc symmetric under a -o- b can be uniquely ex- 
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panded as 



a abc = ic {tbcdOiad + tacdOLbd) + ^Jabc, (25) 



where 



OLda 



1 s 

(Xda 



^Tr[a]S ad + -^-"fdbctbc 



(26) 



(here S ad is the Kronecker delta) and 



labc 



1 



'cab 



°bca) 



3^ [ abc 
1 / 

n \ fa be + Icab + 76 ca J • 



(27) 



Replacing Eq. (19) with Eq. (25) removes the above- 
mentioned discrepancy, because the totally symmetric 
tensor j abc has only ten independent quantities, com- 
pared to 18 in t abc . As for the tensor a, it reduces in 
the static limit to the traceless part of the magnetoelec- 
tric tensor a. But while a becomes origin dependent 
at finite frequencies, 7 a remains origin independent by 
admixing some quadrupolar character. It seems appro- 
priate to interpret the renormalized property tensor a as 
the traceless optical magnetoelectric tensor, and 7 as the 
purely quadrupolar part of cr^ bc . 

We now turn briefly to (T^ bc . A third-rank tensor an- 
tisymmetric in two indices has nine independent compo- 
nents, however, there are 18 quantities on the right-hand 
side of Eq. (20). We therefore replace it with 



a abc 



= ic (ebcdPad - CacdPbd^j , 



where 



Pab = —e bcd (2cr acd - a cda ) 



(28) 



(29) 



— Pab ~\~ ^ C6cd(2^ocd (,cda)- 



Hence natural optical activity, just like gyrotropic bire- 
fringence, is governed by an origin-independent combina- 
tion of magnetoelectric (/?) and quadrupolar (£) terms. 25 
Alternatively, f3 can be interpreted as a renormalized 
magnetoelectriclike tensor, in the same way as a. 

Equations (25) and (28) for a^ bc and a^bc correspond 
to Eqs. (21) and (30) of Ref. 4 while Eqs. (26), (27), 
and (29) express the translationally invariant property 
tensors 5, /3, and 7 as combinations of origin-dependent 
multipole polarizabilities. 



III. EVALUATION OF THE CONDUCTIVITY 

In this section we derive, working in the independent- 
particle approximation, a quantum-mechanical expres- 
sion for o~ abc (uj). The expression, valid for band insula- 
tors, is conveniently written as a sum of two terms, which 



we shall denote by the superscripts (m) and (e). They 
arise, respectively, from the q dependence of the transi- 
tion matrix elements and of the transition energies. 8 

At nonabsorbing frequencies <r abc (uj) is an anti- 
hermitian tensor in the first two indices. The imaginary 
(symmetric) part is given, at T = 0, by the sum of 



Im4 m a l» = ^/[*]E 



(30) 



xlm (Ai nib B n i tac + Ai n _ a B n i_ bc ) 



and 



Imag b » = ir l[dk] 



2e 2 



fi 2 
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E 



,3 

J ln 



n,l v"In (31) 

xd c {Ei+E n )Re(A nl 

,aAl n h) > 

and the real (antisymmetric) part is the sum of 



2e 2 



[dk] 



o.e 



n,l "In-- (32) 
xRc (Ai ntb B n i iac — Ai n _ a B n i^ bc ) 



and 



(33) 



xd c (E t + E n )lm(A nl , a A 



In.b) 



In these expressions the indices n and I run over occu- 
pied (o) and empty (e) bands, respectively, [dk] stands 
for d 3 k/(2n) 3 , d c = d/dk c , and fnoi n — E[—E n . All quan- 
tities in the integrands are labeled by the index k, which 
has been omitted for brevity. The matrix A n i, a = A* ln a , 
known as the Berry connection, is defined as 



A n l,a = i(u n \d a Ul) 



(34) 



and the matrix B n i, ac = —B\ n ac has both orbital and 



spin contributions, 



B nlac = B {o " h) +B\ . 

ni,ac nl,ac ' nl,ac ' 



j(spin) 



(35) 



jiven by 

B { n °£l = Yh K u n\(9 a H)\d cUl ) (d c u n \(d a H)\ Ul )} (36) 



and 



B 



(spin) 
nl.ac 



e a bc(u n \Sb\ui) 



(37) 



where u„k is a cell-periodic Bloch state, is related to 
the crystal Hamiltonian % by e~* kr 'He lk ' r , and m e is the 
electron mass. 

The energy (e) terms have purely orbital character, 
while the matrix element (m) terms have both orbital 
and spin components. It can be verified that the spin 
part of Eq. (30) does not contribute to Eq. (27), consis- 
tent with the fact that "f a bc is a purely orbital (electric- 
quadrupolar) quantity. 
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A. Derivation 

The derivation of the equations given above proceeds 
as follows. We first evaluate the absorptive (Hcrmitian) 
part of a abc , and then insert its symmetric and antisym- 
metric parts into the Kramers-Kronig relations 



Im<7 a6c (w ) 



30 Rea abc (uj) 



duj (38) 



and 



Recr a6c (w ) = -P / duj, (39) 

J-oo OJ-OJ 

respectively. 

The Kubo-Greenwood formula for the absorptive part 
of the conductivity at finite ui and q reads 

^ab(q^) = ^- I [dk] (/n,k-q/2 - flM+q/2) 

nl 

X (VVk-q/2|4 (q) l^,k+q/2) (V'/,k+q/2 l^b(q) |-0n,k-q/2) 

x 5[u - w;„ k (q)] , 

(40) 

where /„k± q /2 is the occupation factor of the Bloch state 
ipnk±q/2 with eigenenergy £„ k±q / 2 , 

^nk(q) = -B;,k+q/2 - ^n,k-q/2, (41) 

and I(q) is related to the velocity and spin operators by 
I(q) = ± + — (S x q) e ^ r . (42) 

Z 777 g 

Equation (40) reduces in the limit q — > to the familiar 
expression for the optical conductivity in the electric- 
dipole approximation. 26 It can be derived starting from 
the interaction Hamiltonian 

ffi = -^(A-v + vA) + — (VxA)-S. (43) 
2c m e c 

Up to terms linear in q, the optical matrix element 
(V'n,k-q/2|4(q)IV';,k+q/2) may be replaced by 



#nik,a(q) = (W„,k-q/2|fa(k) (S X q) a Kk+q/2) 

VClp 



B nLa + B nl,a C q C 



(44) 



where v(k) = e _lk ' r ve jk ' r . Using the relation 13 ftv(k) = 
d a H k together with Eq. (34), the expansion coefficients 
in the second line are found to be 

B nl,a = ( u n\v a \ui) = iu n iA nlta + ^5i n d a Ei (45) 

and Eqs. (35)-(37) for B n ^ ac . 



We are now ready to calculate <J^ bc by differentiating 
Eq. (40) with respect to q c . Because we assume an insu- 
lator at T = 0, 27 the derivative acts only on the transi- 
tion matrix elements and on the S function selecting the 
transition energies, not on the occupation factors. Using 
Eq. (44) for the matrix elements [note that the second, 
intraband, term in Eq. (45) does not contribute in insu- 
lators], together with 



d_ 



5 [u - w*„k(q)] 



q=0 



(46) 



- — S'(w - uj lnk (0))d c (E lk + E nk ) 



and inserting the result for the symmetric and antisym- 
metric parts of a^ bc into Eq. (38) and Eq. (39) respec- 
tively, one easily obtains Eqs. (30)-(33). 



B. Static limit 

In the limit w — > the ME tensors xTb an d Xba become 
identical, and as a result a^ bc [Eq. (20)] vanishes. As for 
a abc we noted in Sec. II C that its dc limit is governed 
by 5(0), the traceless part of the static ME polarizability 
tensor a(0). Since our calculation of crf 6c only included 
the purely electronic response to the optical fields, we 
should recover in that limit the frozen- ion part of 5(0). 

We will focus here on the orbital contribution to crj bc , 
and compare it with the band-theory expression obtained 
in Refs. 20 and 21 for the frozen-ion orbital ME tensor. 
The corresponding proof for the spin contribution is ele- 
mentary. 

We begin by recasting the orbital part of Eqs. (30) and 
(31) at u! = in a form where empty states do not appear 
explicitly This is done in Appendix B, where we obtain 

2 r ° 

Im 4°abc(°) = \\ [ dk ] ^2~Rc{{d a u n \d c u m ){u m \d b u n ) 

nm 

+ {d b u n \d c u m ){u m \d a u n )^ 

/o 
[dk] { [lm(d c u n \d a (H + E n )\d £b u n ) -a*+c] 
n 

+ b <-> a j, 



+ 



(47) 



where the covariant field derivative \ds b u n ) is given by 
Eq. (B5). Equation (47) can now be compared with 
Eq. (C.2) of Ref. 20 for the static ME tensor, which reads 



/o 
[dk] ^2 Im{dbU n \d c {H + E n )\d £d u n ). 



(48) 
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FIG. 1. (Color online) The xxy component of the gyrotropic 
birefringence tensor Imcrf i)c , and the xyz component of the 
natural optical activity tensor Recr^, calculated for the 
tight-binding model described in the text as a function of fre- 
quency. Solid lines: extrapolation from calculations on finite 
crystallites. Dashed lines: calculations on periodic crystals 
using the fc-space formulas derived in this work. The verti- 
cal dotted line indicates the frequency corresponding to the 
direct band gap. 
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FIG. 2. (Color online) The xxy component of Ima^ bc (uj), 
calculated for the tight-binding model described in the text 
as a function of the parameter <p. Solid lines: extrapolation 
from calculations on finite crystallites using Eq. (49). Dashed 
lines: calculations on periodic crystals using Eqs. (30) and 
(31). Dotted lines: same as the dashed lines, but ommiting 
the contribution coming from Eq. (31). 



It is easily verified that inserting Eq. (48) into Eq. (19) 
at to — yields Eq. (47), which proves the result. 

IV. NUMERICAL RESULTS 

In order to check the expressions derived in the pre- 
vious section, we have carried out numerical tests com- 
paring calculations done under periodic boundary condi- 
tions against reference calculations on finite crystallites. 
We chose for our tests the tight-binding model of Ref. 20. 
This is a spinless model on a 2 x 2 x 2 cubic lattice, where 
V symmetry is broken by assigning random on-site ener- 
gies and T symmetry is broken by complex first-neighbor 
hoppings. The model parameters in Table A.l of Ref. 20 
were used (one of the complex hopping phases, labeled 
ip therein, shall be used as a control parameter), and the 
two lowest bands were treated as occupied. 

The tensor components Im erf and Re a^ yz were eval- 
uated at nonabsorbing frequencies. The calculations on 
periodic samples were done on a 30 x 30 x 30 mesh of 
k points using Eqs. (30)-(36), together with the sum- 
over-states formula for Vk|itnk)- 20 For the calculations 
on finite samples we used Eqs. (19) and (20), 

Im erf xy = 2cRea xz + ulmjxxy = 2ca xz - iuj xxy (49) 

and 

Re a xyz = -clm(/3 xx + p yy ) + u)Re£ xyz 

— ic^ftxx 0yy) ^£,xyz-, 

together with Eqs. (21)-(24) and (Al)-(A4) for the mag- 
netoelectric (a, (3) and quadrupolar (7, £) tensors. We 



chose cubic samples containing L x L x L unit cells, with 
L = 1, 2, 3, 4, and then extrapolated the calculated val- 
ues to L — > 00. 20 

Figure 1 shows as solid (dashed) lines the frequency 
dependence of Im <J xxy and Re <J^ yz for finite (periodic) 
samples, with the parameter <p set to it. The natural opti- 
cal activity spectrum starts off at zero and increases with 
frequency, exhibiting a resonant behavior as the mini- 
mum direct gap, denoted by the vertical dashed line, is 
approached. The ME optical spectrum displays a similar 
behavior, except that it remains finite as tu goes to zero. 
The excellent agreement between solid and dashed lines 
demonstrates the correctness of the fc-space formulas. 

Next we discuss a number of additional numerical 
tests where we investigate in more detail the behavior 
of lma xxy . In these tests the frequency was kept fixed, 
and the parameter <p was scanned over the range [0, 2tt]. 

In Figure 2 we plot Im erf versus <p for two frequen- 
cies, lu — and fkj = 1. As before, solid and dashed lines 
represent calculations on finite and periodic samples re- 
spectively. In addition, we show as dotted lines the result 
of a periodic-sample calculation using only the matrix el- 
ement (m) term, Eq. (30), i.e., omitting the energy (e) 
term, Eq. (31). We see that the energy term gives a small 
but visible contribution, which must be included in order 
to find agreement with the finite-sample calculation. 

We now turn to the decomposition of Im erf^ accord- 
ing to Eq. (49), into magnetoelectric and quadrupolar 
parts. They are plotted separately in Fig. 3 for hw = 1 
and L — 4. We chose a specific L because a and uij 
are origin-dependent quantities, and it is therefore not 
meaningful to extrapolate them separately to L — > 00. 
The dashed lines show how each of them changes when 
the position of the sample is shifted. The change in a zz is 
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FIG. 3. Origin-dependence of the bare magnetoelectric (up- 
per panel) and quadrupolar (lower panel) polarizabilities ap- 
pearing on the right-hand side of Eq. (49), calculated at 
hid = 1 for a finite sample (L = 4) of the model used in 
Fig. 2. Solid lines: the center of the sample is placed at the 
origin. Dashed lines: the sample is displaced by r = (1, 1, 1), 
in units of the lattice constant of the 2x2x2 cubic cell. 



exactly compensated by the change in uj xxy , so that the 
resulting Im o xxy remains the same to machine precision, 
demonstrating its translational invariance. 

An alternative decomposition of Imcr^ is given by 



Eq. (25): 



(51) 



Unlike the bare property tensors a and uij appearing 
in Eq. (49), the renormalized magnetoelectriclike and 
purely quadrupolar tensors a and urj are origin inde- 
pendent and hence separately well defined for periodic 
samples. Figure 4 shows as dashed (solid) lines their val- 
ues calculated for periodic (finite) samples from the first 
(second) equality in Eqs. (26) and (27). Because a re- 
duces to the traceless part of a as lj — > 0, we can directly 
compare the curve for a xz (0) with a fc-space calculation 
of a xz (0) using the formula derived in Refs. 20 and 21 
(open circles). The precise agreement confirms numeri- 
cally the analysis of Sec. IIIB. 



FIG. 4. (Color online) Translationally invariant decomposi- 
tion [Eq. (51)] of the curves in Fig. 2 into magnetoelectriclike 
(upper panel) and purely quadrupolar (lower panel) contri- 
butions. Solid lines: extrapolation from calculations on finite 
crystallites. Dashed lines: fc-space calculations on periodic 
crystals. In the static limit the tensor 5 reduces to the trace- 
less part of the magnetoelectric polarizability a, and the open 
circles show a xz (0) calculated in k space according to Refs. 20 
and 21. 



V. SUMMARY AND OUTLOOK 

In this work we investigated spatial-dispersion optical 
effects in insulators. The main result is a band-theory ex- 
pression for a a b c (u)), the spatially dispersive optical con- 
ductivity. Special attention was given to the T-odd part 
of this tensor, which is nonzero in magnetoelectric crys- 
tals, and comprises magnetoelectriclike (a a b) and purely 
quadrupolar ("f a bc) contributions. We showed that each 
of them consists of a translationally invariant combina- 
tion of separately origin dependent molecular polarizabil- 
ity tensors. 

The magnetoelectriclike tensor a a b has both spin 
and orbital contributions, and the expression for the 
orbital part generalizes to finite frequencies the re- 
cently developed band theory of orbital magnetoelectric 



response 



20,21 



The generalization is, however, not com- 
plete, as the tensor a a b(u) is traceless, and therefore does 
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not include the isotropic ME coupling. The reason why 
the latter is not recovered from the present formalism is 
that our starting point is the current response of an infi- 
nite medium to an electromagnetic wave while the trace 
of the ME tensor, known as the axion contribution, only 
affects electrodynamics at boundaries. 4 ' 21 The calcula- 
tion of the axion piece at finite frequencies remains an 
open problem. 

The bulk expression for cr a h c (w) at transparent frequen- 
cies was validated by performing numerical calculations 
on a tight-binding model, and comparing against refer- 
ence calculations done on finite samples. The quantities 
needed to evaluate that expression are the occupied and 
empty energy eigenvalues and their fc-space gradients, the 
off-diagonal Berry connection matrix Eq. (34), and the 
orbital and spin matrices Eqs. (36) and (37). The evalua- 
tion of all these objects in a first-principles context can be 
done efficiently by mapping the electronic structure onto 
localized Wannier orbitals, and then using the technique 
of Wannier interpolation. 28 This approach has already 
been used to compute the magnetic circular dichroism 
spectrum of ferromagnets. 29 

First-principles calculations of the optical spectrum of 
solids beyond the electric-dipole approximation are still 
in their infancy. We hope that the formalism introduced 
in this work will be useful for carrying out realistic cal- 
culations of spatial-dispersion phenomena in the optical 
range, including natural optical activity, gyrotropic bire- 
fringence, and directional dichroism. 
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ImxX = ^-Im [<n|r |i)<i|(rxv) 6 |n>] = -Im X , 

n.l 



mc 
ba 



and the quadrupolar polarizability reads 



(A2) 



n.l 



Rex^ c = E ^ Re [{n\r a \l){l\nr c \n)] = ^X% ba , 

(A3) 



n.l 



lm Xabc = E ~ Im [(n\r a \l)(l\r b r c \n)] = -Im^. 

(A4) 



Appendix B: Derivation of Eq. (47) 

In order to derive Eq. (47), we drop the spin contribu- 
tion, Eq. (37), from Eqs. (30) and (31) and rewrite the 
orbital contribution at u = as 

2 r- 

Im4°abc(0) = \ J [<*] KCabc + D abc ) + (C bac + D bac )\ , 

(Bl) 

where 



c abc = j2 Ri 



(ui\d b u n ) 



Ei — E„ 



[u„\(d a H)\d c ui)- 



(B2) 



{d c u n \(d a H)\ui) 



and 



D abc = E d ° f^T^ ^ Re{(d a u n \ui)(ui\d b u n )}. (B3) 



Appendix A: Quantum-mechanical expressions for 
the polarizability tensors 

In this appendix we list the quantum-mechanical ex- 
pressions for the frequency^dependent polarizability ten- 
sors x cm , x mc , x q , and x q 01 bounded samples. They 
have been used to produce the reference results (solid 
lines) in Figs. 2-4. 

We provide the single-particle version of the formulas 
in the lossless regime, which is the form used in Sec. IV. 
A many-body derivation can be found in Rcf. 7, and the 
modifications needed to describe absorption are discussed 
in Refs. 6 and 7. 

Defining Zi n = (Vh/2e 2 )(uf n — w 2 ), where V is the 
system volume, the orbital contribution to the magneto- 
electric tensor reads 



RexX = Y c E |; Re MralOWr x v) 6 |n>] = Re X £f 

n,l 

(Al) 



We will use repeatedly the identity 20 

d c {H-E l )\u l ) = {E l -H)\d c ui), (B4) 

as well as the following expression for the field derivative 
of a valence-band state projected onto the conduction 
bands 20 

\d £b u n ) = - l eY j ^^\d b u n ). (B5) 

We start by using Eq. (B4) to eliminate d c Ei from 
Eq. (B3), 



Dabc 



u n \{d c H)\ui)(ui\d b u n ) 



Ei - E„ 



+ 



(d a u n \H - Ei\d c ui){ui\d b u n ) 



d c E n 

El - E n 



Ei - E n 
{d a u n \ui){ui\d b u n 



(B6) 
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and then use Eq. (B5) twice to find 
1 ° 

Dabc = — ^2lm(d a u n \d c (H + E n )\d £b u n ) 



n,l 



Ei\d c ui)(ui\d b u. 



Ei - E n 



(B7) 



Now write H - Ei as (H - E n ) + (E n - E{) and use 
Eq. (B4), 

1 ° 

Dabc = T abc - - ^2lm(d a u n \d c (H + E n )\d £b u n ) 



+ E Re 



n,l 



(u n \d a (E n - H)\d c ui)(ui\d b u n ) 

El — E n 



(B8) 



where we defined 



Tabc = -^2Rc(d a u n \d c ui)(ui\d b u n ). (B9) 

n,l 

One term in Eq. (B8) exactly cancels the first term 
in Eq. (B2). For the remainder we use (u n \d c ui) = 



— (d c u n \ui) once and then Eq. (B5) twice, yielding 

C a bc ~t~ D a bc — T abc 

1 ° ~ 

+ -^Im{(<9 c u„|<9 Q (iJ + E n )\d £b u n ) - a <H- c| . 

n 

(BIO) 

In order to eliminate the sum over empty states in T abc 
we need to combine C a bc + D abc with C bac + D bac , as in 
Eq. (Bl). We therefore consider 

o.e 

T a bc + T bac = - ^2Rc{{d a u n \d c ui}{ui\d b u n } + a o &} 

n,/ 
o 

= - ^2~R.c{(d a u n \(d c Q)\d b u n )} 

n 

o 

u n \d c u m ){u m \d b u n ) 

nm 

+ (d b u n \d c u m ){u m \d a u n )Y 

(Bll) 

where Q = |«j)<«z I = 1 - Em l"m)(«m|- Collecting 
terms, we arrive at Eq. (47). 
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